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ABSTRACT 


Piwi-interacting RNA (piRNA) plays an important role 
in the gonadal development and maintenance of 
Teleostei. In this study, piRNA libraries derived from 
the adult gonads of Japanese flounder (Paralichthys 
olivaceus) were generated using next-generation 
sequencing technology. Using zebrafish piRNAs as 
a reference, 5 865 unique candidate piRNAs were 
identified; 289 candidate piRNA clusters (PRCs) 
were generated from the above piRNAs. Among the 
isolated candidate PRCs, a total of 38 ovary-specific, 
45 ovary-bias, 24 testis-specific, and 131 testis-bias 
PRCs were found. The relative expression levels of 
seven PRCs were validated through quantitative 
reverse transcription-polymerase chain reaction. The 
results of this study will help facilitate exploration of 
the development and maintenance of the phenotypic 
sex mechanism in P. olivaceus. 
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INTRODUCTION 


In mammals, the corresponding protein of piwi-interacting RNA 
(piIRNA) interacts with a class of germ cell-restricted small 
RNAs (26-31 nucleotides in length) in the gonads. piRNAs 
show distinctive localization patterns in the genome and are 
predominantly grouped into 20—90-kilobase clusters, whereas 
the long stretches of small RNAs are derived from only one 
strand. In humans and rats, with major clusters occurring in 
syntenic locations, both the abundance of piRNAs in germline 
cells and male sterility of Piwi mutants suggest the possible role 
of piRNAs in gametogenesis (Aravin et al., 2006; Girard et al., 
2006). In zebrafish (Danio rerio), the piwi gene is expressed in 
both the testis and ovary and is also a component of a germline 
structure called “nuage”. The loss of piwi function results in a 
progressive loss of germ cells due to apoptosis during larval 
development. Many of these small RNAs are derived from 
transposons; therefore, piRNAs are possibly involved in the 
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silencing of repetitive elements in gonads (Houwing et al., 
2007). In general, genome-wide piRNA expression and function 
have been examined in fruit flies (Drosophilidae) (Aravin et al., 
2007), humans (Homo sapiens) (Sigurdsson et al., 2012), mice 
(Mus musculus) (Zheng et al., 2010), zebrafish (Houwing et al., 
2007; Huang et al., 2011), and planaria (Prosorhochmus 
claparedii) (Palakodeti et al., 2008). Furthermore, piRNAs are 
found in clusters throughout the genome, and each cluster may 
contain as few as 10 or up to thousands of piRNAs. The 
clustering of piRNAs is highly conserved across species, 
whereas the sequences are not (Malone et al., 2009). 

The gender of Japanese flounder (Paralichthys olivaceus) is 
easily altered by water temperature during the sex 
determination stage (Wen et al., 2014). Female P. olivaceus 
grow faster than males; thus, investigations on the mechanisms 
of sex differentiation can be potentially useful in increasing 
production of farmed P. olivaceus. Study on sex manipulation in 
P. olivaceus has been a topic of research for some time. Sex- 
related genes have been studied to elucidate the sexual 
molecular mechanism in flounder, and include the cyp19a gene, 
which was first cloned and found to be predominantly 
expressed in female flounders (Kitano et al., 1999); some 
transcription factors, such as foxl2 (Yamaguchi et al., 2007), 
dmrt1 (Wen et al., 2014), and dmrt4 (Wen et al., 2009), and 
other sex-related genes, such as cyp17 (Ma et al., 2012). 
However, due to the limited data on these genes, information 
regarding their expression and regulation mechanisms is also 
incomplete. As a representative of non-coding RNA, 
approximately 20.0% and 13.1% of isolated microRNAs 
(miRNAs) are preferentially expressed in the testis and ovary, 
respectively, which indicates that miRNAs might play a key role 
in the regulation of gene expression in P. olivaceus gonads (Gu 
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et al., 2014). However, the differential expression and function 
of piRNAs in fish sex development and maintenance remain 
poorly investigated. Identifying the expression patterns of 
piRNAs in P. olivaceus gonads is the first step in clarifying their 
function in gonadal development. Thus, this study aimed to 
identify candidate piRNAs in sexually mature gonads of P 
olivaceus through next-generation sequencing, which provides 
not only sequences of low-abundance piRNAs, but also 
quantitative data on the frequency of sequencing reads that 
correspond to the abundance of piRNAs in the ovary and testis. 


MATERIALS AND METHODS 


Fish maintenance 

One-year-old P. olivaceus fish with average individual body 
weights of 450-500 g were obtained from the Beidaihe Central 
Aquaculture Experiment Station, Chinese Academy of Fishery 
Sciences. All fish were acclimatized to indoor recirculation 
systems at the Shanghai Ocean University for two weeks. The 
gonads were examined through histological sectioning before 
RNA extraction, and only healthy fish were chosen for RNA 
isolation. All experiments were performed according to the 
Experimental Animal Management Law of China and as 
approved by the Animal Ethics Committee of Shanghai Ocean 
University. 


Histology and RNA isolation 

All fish were killed by overdose with tricaine methanesulfonate 
(MS-222, 200 mg/L). Mature ovaries and testes were collected 
separately from five females and five males. For histology, the 
gonads of females and males were fixed with Bouin’s solution, 
then embedded in paraffin (Sigma Aldrich, http: 
/Iwww.sigmaaldrich.com). Cross sections were cut at 5 ym and 
stained with hematoxylin and eosin. Tissues were isolated, 
immediately frozen in liquid nitrogen, and stored at —80°C. Total 
RNA was extracted using a mirVana™ miRNA isolation kit 
(Ambion; http: //www.thermofisher.com/cn/zh/home/brands/ 
invitrogen/ambion.html) in accordance with the manufacturer’s 
instructions. Each testis or ovary was separately subjected to 
RNA extraction. Five RNA samples (1.5 ug per sample) derived 
from the same type of gonad were pooled for library 
construction and subsequent sequencing. The quantity and 
purity of total RNA were monitored using Nanodrop ND-1000 
(Nanodrop Technologies; http: /Awww.thermoscientific.com/content/ 
tfs/en/product/nanodrop-lite-spectrophotometer.html) and denaturing 
gel electrophoresis. The total RNA was stored at —80°C. 


Small RNA library construction and sequencing 

Two small complementary DNA (cDNA) libraries were 
generated from the mixture of total RNAs of mature testes 
and ovaries, respectively. A combination of RNA samples 
was used to avoid possible expression bias of piRNAs in 
different individuals. A small RNA fraction (15-31 nt) was 
isolated using 15% polyacrylamide/8 mol/L urea/0.5*Tris/borate/ 
ethylenediaminetetraacetic acid gel electrophoresis and ligated 
with proprietary 5' (GTTCAGAGTTCTACAGTCCGACGATC) 
and 3' (TCGTATGCCGTCTTCTGCTTGT) adaptors using T4 
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RNA ligase. The short RNAs were then converted to cDNA 
through reverse transcription polymerase chain reaction (RT- 
PCR). The results were verified through Pico green staining and 
fluorospectrophotometry and then quantified using Agilent 2100; 
afterward, the cDNA libraries were sequenced using HiSeq2000 
(Illumina; http: //www.illumina.com/) in accordance with the 
manufacturer's recommended protocols for small RNA 
sequencing (NEXTflex'” Small RNA-seq kits; BIOO Scientific 
Corp; http: /www.biooscientific.com/). 


Bioinformatics analysis of sequencing data 

Raw sequences were filtered to remove low-quality reads 
(sequences with more than 10% sequencing errors), and those 
that passed the quality filter were trimmed to remove the 
adaptor sequences. Sequences with 15-31 nt in the two 
combined libraries were then selected as raw small RNA 
sequences. The raw small RNA sequences were mapped to the 
piRNA database (Lakshmi & Agrawal, 2008) with zebrafish 
piRNAs used as the templates. The 5 865 candidate piRNAs 
are listed in Supplemental Table 1. Highly conserved unique 
piRNAs with different lengths were placed in the same piRNA 
clusters (PRCs), and total reads of all unique piRNAs were 
considered as the reads of the PRCs. The expression levels of 
the PRCs were estimated using RPKM (reads per kilobase of 
exon model per million mapped reads). Differential expression 
analyses of the two samples were performed using the DESeq 
R package (Anders et al., 2013). P-values were adjusted using 
the Benjamini and Hochberg method. A corrected P-value of 
0.05 was set as the threshold for significantly differential 
expressions. The P-values were also adjusted using Q-values; 
Q-value <0.01 and logs (fold change) >1 were set as the 
threshold for significantly differential expressions by default 
(Storey, 2003). 


Table 1 Primers used for piRNA qRT-PCR analysis 


Primer name Sequence(5’-3’) 





piR_0048450 GATGAACCGAACGCCGGGTTAA 
piR_0006070 AAGATGGTGAACTATGCCTGGG 
piR_0019030 GCAAATCGGTCGTCCGACCTGG 
piR_0015958 GATGAACCGAACGCCGGGTTAA 
piR_0055443 TATCATCCAGGTACACAAAGACA 
piR_0056629 TCCATTTTCAGGGCTAGTTGATT 
piR_0050867 TGATGATAGGTGTGAGTGCTACA 
Po5Ss CCATACCACCCTGAACAC 

Po5Sa CGGTCTCCCATCCAAGTA 


Quantitative RT-PCR analysis of piRNAs 

The expression profiles of nine PRCs were investigated with 
real time RT-PCR using a miScript reverse transcription kit and 
miScript SYBR Green PCR kit (Qiagen: http: //www. 
transmedchina.com/) in accordance with the manufacturers’ 
instructions. Quantitative RT-PCR (qRT-PCR) was analyzed 
using iCycler (Bio-Rad; http: //www.bio-rad.com/). The 20 uL 
reaction mixture consisted of 10 uL of SYBR Premix Taq (2x), 


0.4 uL of piRNA-specific forward primer (10 mol/L), 0.4 uL of 
miScript universal primer (10 pmol/L), and 1 uL of PCR 
template (CDNA). The amplification protocol was as follows: 
initial denaturation and enzyme activation for 3 min at 95°C, 
followed by 40 cycles for 5 sec at 95°C, and for 20 sec at 60°C. 
Fluorescence was read at the end of each round of 
amplification. Each assay was performed in duplicate. The 
relative expression levels of one piRNA in the ovary and testis 
were determined using the (1+E7)°/(1+ER)“™ method, and 
different piRNAs in the same tissue were determined using the 
2°) method (Pfaffl, 2001). The 5S small ribosomal RNA of P. 
olivaceus was used as the internal control, and 
piRNA_0002790 was used as the calibrator. The primers of 
each piRNA are listed in Table 1, with the most conserved 


sequence of each PRC selected as the primer. The groups 
were compared via one-way analysis of variance (ANOVA) and 
Tukey's test to identify statistically distinct groups. Significant 
differences were accepted using P<0.05. 


RESULTS 


Histological observation and identification of piRNAs of gonads 
In the ovary of P. olivaceus, only mature oocytes or those close 
to maturity were revealed. The previtellogenic oocytes 
displayed strongly basophilic ooplasms, while mature oocytes 
were strongly eosinophilic. In the testis, tubules were filled with 
spermatozoids, while the numbers of spermatogonia and 
spermatocytes were reduced (Figure 1). 
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Figure 1 Histological features of the ovary and testis of one-year-old P. olivaceus 
A: Ovary. 1, Stage Ill; 2, Stage VI; 3, Stage V. B: Testis. b, blood vessel; f: flagella; arrow indicates interstitial cells. 


Combined RNA samples from different ovary and testis 
samples were used for sequencing. A total of 24 399 395 and 
26 417 330 raw reads were obtained from the ovary and testis 
libraries, respectively. After filtering low-quality sequences, 
empty adaptors, and single-read sequences, 12 246 743 
(50.19%) and 13 260 806 (50.19%) clean reads were used for 
further sequence analyses. After assembly and BLAST 
searching in the piRNA database, 5 181 and 3 915 unique 
piRNA sequences showing high homology with zebrafish were 
obtained from the two libraries. All unique piRNAs were 
grouped into 289 piRNA clusters (PRCs), where all unique 
piRNAs in the same PRC shared a conserved sequence with 
varying lengths. The remaining nucleotides were similar to other 
types of RNAs, including ncRNA, tRNA, miRNA, rRNA, snRNA, 
and snoRNA. The size distribution of small RNAs for 
sequencing was similar between the two libraries. The majority 
of these small RNAs increased from 18 nt to 30 nt, and lengths 
ranged from 23 nt to 28 nt (Figure 2A). 


Expression patterns of PRCs in P. olivaceus gonads 

The relative abundance of various PRCs was determined using 
the small RNA deep sequencing approach by calculating 
sequencing frequency. In the ovary and testis libraries, the 
identified PRCs exhibited a broad range of expression levels. 
The most abundant PRC in the mature gonads was 
piR_0000766 (which displayed more than 12 000 reads in the 
ovary and testis), followed by piR_0034879 (10 766 reads). In 


contrast, 87 piRNAs displayed less than 10 reads. These 
results indicate that expression varied significantly among 
different piRNA families. 


Differential expression of PRCs in the ovary and testis 

On the basis of changes in the relative piRNA abundance 
between the two libraries, we found that logs'*"*s"s°""™) of the 
expressed PRCs varied from -7.65 to 4.89. Approximately 
45.32% and 15.57% of isolated PRCs were preferentially 
expressed in the testis and ovary, respectively. Results showed 
a total of 289 PRCs, including 38 ovary-specific, 45 ovary-bias, 
24 testis-specific and 131 testis-bias PRCs (Figure 2B). Seven 
PRCs (piR_0048450, piR_0006070, piR_0019030, 
piR_0015958, piR_0055443, piR_0056629, and piR_0050867) 
with different reads obtained from sequencing results were 
selected and quantified using real-time RT-PCR to validate the 
results obtained from the Illumina deep sequencing data. These 
PRCs were also homologues, according to the piBASE 
database (http: //www.regulatoryrna.org/database/piRNA/). 
Three individual experiments using different RNA pools were 
performed to generate reliable and repeatable results. Results 
showed that PRCs, including piRNR_0048450, piRNA_0006070, 
piRNA_0019030, and piRNA_0015958, exhibited high 
expression levels in both the testis and ovary, though were 
more highly expressed in the testis (P<0.05). Three other PRCs 
(pIRNA_0055443, piRNA_0056629, and piRNA_0050867) 
were primarily detected in the ovary (P<0.05). These results are 
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consistent with the sequencing data (Figure 3). 
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Figure 2 Length range and differential expression of piRNA in 
the gonads of P. olivaceus 

A: Length distribution of sequenced small RNAs from the ovary and 
testis of P. olivaceus. X-axis denotes the piRNA length of 26-31 
nucleotides; Y-axis represents the percentage of piRNA. B: Scatter plot 
expression level of piRNA in testis versus ovary. Genes are listed in 
Supplemental Table 2. 


DISCUSSION 


As the largest group of vertebrates, teleosts exhibit remarkably 
varied sexuality (Godwin, 2010). Genetic, environmental, and 
social factors can affect sexual development in a large number 
of fish species at different developmental stages (Heule et al., 
2014). Paralichthys olivaceus, which is a common flatfish raised 
in China, Korea, and Japan, exhibits a temperature-sensitive 
mechanism of gender determination and differentiation. Artificial 
breeding and selection is performed to facilitate the production 
of P. olivaceus, and sex control is regarded as an important 
technique for fish production (Zhu et al., 2006). Although the 
sexual development of P. olivaceus needs to be properly 
understood, information concerning the involved mechanism 
remains insufficient and unclear. One-year-old Japanese 
flounder are suitable for supply to the market. The gonad is 
mature and easy to collect, observe and study in regards to 
gene expression. 

Three processes take place during sexual development, that 
is, sex determination, differentiation, and maintenance. The first 
process determines whether the bipotential primordium will 
develop into a testis or ovary. The second process involves the 
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Figure 3 Relationship between relative expression levels of 
piRNA validated by qRT-PCR and fold changes derived using 

small RNA sequencing 

Expression levels of piRNA measured with qRT-PCR are depicted with 
column charts, and triangles or circles represent the fold changes in 
sequencing. Error bars represent one standard deviation of three 
different biological replicates. Statistical significance is reported for each 
piRNA. F: female; M: male; *: P<0.05. 


actual development of the testes or ovaries from the 
undifferentiated gonads. The final process ensures the 
reproductive function of the ovary and testis (Freedberg & 
Taylor, 2007; Hayes, 1998; Spencer & Janzen, 2014). The 
protein-coding genes involved in the regulatory network of 
sexual development described so far originate from mammalian 
or Teleostei studies and can be divided into different functional 
groups. The dmy gene has been identified in medaka and plays 
a key role in gender development (Matsuda et al., 2002); other 
genes, including dmrt3 (Li et al., 2008), dmrt6 (Zhang et al., 
2014), gsdf (Shibata et al., 2010), 118-hydroxylase cyp11b 
(Schiffer et al., 2015), and 5a-reductase 1, 2, and 3 (Aquila et 
al., 2015), are also involved in sexual development. In addition, 
some genes, such as dmy/dmrt1 in medaka, possess different 
biological functions in different strains (Matsuda et al., 2002). 
Plasticity is the main characteristic of sexual development, 
especially for fish. Studies on noncoding genes, including 
piRNA, microRNA, and IncRNA, can provide insights into 
sexual identification in fish (Mishima et al., 2008; Gu et al., 2014; 
Jing et al., 2014; Xiao et al., 2014). Gender-specific 
organization of piRNA has been discovered in Drosophila (Li et 
al., 2009; Malone et al., 2009), as have differences between 
male and female piRNAs in zebrafish (Houwing et al., 2007). 
Through microfluidic chip analysis, piRNAs have been found to 
exhibit quantitative differences, ranging from 1.58 to 64 019.26. 
A subset comprising more than 20% of piRNA exhibited evident 
differential expression between the testis and ovary (Houwing 
et al., 2007). We also observed gender-based differential and 
specific organization of piRNA in P. olivaceus. Results suggest 


that the gender-related organization of piRNA might be 
conserved in vertebrates and flies, and might possess multiple 
functions in organisms. Piwi-antisense piRNA complexes 
mediate cleavage of sense piRNA precursors and transposon 
(or protein-coding) transcripts, which silence transposon and 
gene expression at the post-transcriptional level (Weick & 
Miska, 2014). Therefore, study on piRNA is important for 
illuminating the function of sex development-related protein- 
coding genes. 

piRNAs are involved in the sexual development of many 
vertebrates, particularly mammals and teleosts (Aravin et al., 
2007; Girard et al., 2006; Houwing et al., 2007). However, the 
related mechanism remains only partially understood because 
of the lack of information on piRNA expression patterns in the 
gonads of non-model organisms. Furthermore, identification of 
piRNAs is difficult because they lack conservative secondary 
structure motifs, such as the ones found in miRNA (Zhang et al., 
2010). Thus, we examined all potential piRNAs in both male 
and female gonads through deep sequencing technology, 
followed by BLAST searching the piBANK database with the 
zebrafish piRNA library as a reference. With a high frequency of 
24-28 nt RNA, the deep sequencing results suggested that 
piRNAs are present in the male and female gonads of P. 
olivaceus. These findings are consistent with piRNA length 
distribution in zebrafish. Mammals, including humans, mice, 
and rats, have a major distribution interval of 27—31 nt (http: 
//pirnabank.ibab.ac.in/stats.html). The crystal structure of the 
MID domain is obtained from a Piwi/Argonaute protein, which 
plays a central role in RNA silencing process as essential 
catalytic components of the RNA-induced silencing complex 
(RISC). The RISC complex is responsible for the gene silencing 
phenomenon known as RNA interference. MID docking 
experiments showed its capability to specifically recognize the 
5'-uridine (1U-bias) of piRNAs (Cora et al., 2014). In this study, 
the ratios of A, T (U), C, and G in the first nucleotides of the 5 
865 unique piRNAs of P. olivaceus were 23.24%, 30.75%, 
22.95%, and 23.05%, respectively, which also characterized 
1U-bias of the molecule. The expression levels of PRC in P. 
olivaceus presented a broad range. For example, more than 
several thousand reads of piR_0000766 were found in both 
male and female piRNA libraries, whereas specific infrequent 
PRCs showed several reads only. We suggest that the 
imbalance in the expression levels of different piRNAs might 
correspond to the different abundance levels of their target 
genes. Hence, the relationship between the expression levels 
of PRCs and their target genes should be considered in future 
studies. 
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